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Abstract 

The connection between anomalous scaling of structure functions (intermittency) and 
numerical methods for turbulence simulations is discussed. It is argued that the com- 
putational work for direct numerical simulations (DNS) of fully developed turbulence 
increases as Re 4 , and not as Re 3 expected from Kolmogorov's theory, where Re is 
a large-scale Reynolds number. Various relations for the moments of acceleration 
and velocity derivatives are derived. An infinite set of exact constraints on dynam- 
ically consistent subgrid models for Large Eddy Simulations (LES) is derived from the 
Navier-Stokes equations, and some problems of principle associated with existing LES 
models are highlighted. 



1 



1 Background 



The theory of turbulence and the development of calculation methods for high-Reynolds- 
number flows became an active research topic around the beginning of the twentieth century. 
This effort yielded many important results of general interest in statistical physics. For in- 
stance, Kolmogorov's work [l]-[3] on turbulence theory formulated the scaling ideas for the 
first time, and Kraichnan [4] proposed the mode coupling approach. However, the "turbu- 
lence problem", lacking a small parameter characterizing the strong nonlinear interactions, 
has turned out to be remarkably difficult — and it remains so today. 

The revolutionary realization of Osborne Reynolds that turbulence theory is a subject 
of statistical hydrodynamics rather than classical hydrodynamics, led almost hundred years 
ago to various elegant and useful phenomenological models based on ideas of kinetic the- 
ory (Prandtl [5], Richardson [6], Kolmogorov [3]), which strongly impacted the engineering 
profession. These heuristic semi-empirical models, based on low-order closures of various 
perturbation expansions, had a somewhat limited range of success and needed adjustable 
parameters, often varying from flow to flow. Nevertheless, the role of these models was— 
and still is — so immense that one can hardly imagine processes in mechanical and chemical 
engineering, aerodynamics and meteorology which do not have their input. 

With the advent of powerful computers, the possibility of accurate numerical simulations, 
directly based on the Navier-Stokes equations, became a reality. Since the introduction of 
spectral methods in the end of sixties [7]- [8], direct numerical simulations (DNS) have become 
a new tool to attack the "turbulence problem". A strategic goal of the DNS has been to 
complement expensive and complicated physical experiments, and their dream is to dispense 
with them altogether. 

The computational power required for DNS is estimated on the basis of Kolmogorov's 
phenomenology that describes turbulent fluctuations filling the interval of wavenumbers 
1/L < ^ < where L and r/ K = LRe^* are the integral and dissipation scales, 

respectively, and Re = u Tms L/v is the Reynolds number based on L and the root-mean- 
square velocity u rms . If we assume that the velocity fluctuations on scales r « rjx are 
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highly damped and cannot contribute to the inertial range dynamics, the effective number 
of degrees of freedom [9] is then (L/t]k) 3 = Re 9 / 4 . This is the minimum number of grid 
points required in DNS for a cubic box of linear dimension L. The required number of time 
steps in the computation is usually proportional to the spatial grid points, so the total com- 
putational work increases as Re 3 . This means that a mere doubling of the Reynolds number 
requires almost an order of magnitude increase of computational work. 

The accuracy of numerical methods is traditionally estimated as follows. The dissipation 
contribution to the equation for turbulent kinetic energy is given by 



Q2 U Q2 y d 2 2 

S = -z/u • = -v lim r ^ ri —Ui(x)ui(x + r) = v lim r ^ v - — S 2j0 (r) oc vE*rf 2 ~ 2 , 

where the order of magnitude estimate in the last step comes from Kolmogorov's phenomenol- 
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ogy. For this case, £ 2 = 2/3 and we have r\ K = We then have the familiar estimate 
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i] K S3 LRe~i, mentioned earlier. Thus, to accurately describe the flow, one has to simply 
account for fluctuations on the scales r > t]k by choosing the computational mesh size to be 
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A = ar] K S3 aLi?e~4, (1) 

where a = const = 0(1). On this mesh, the velocity derivative is defined as 

u{x + A) — u{x) = du{x) v 1 d n u{x) ! 

A dx ^ 2 n\ dx n { ' 

Now, in Kolmogorov's turbulence, (d x u) rms = \J (d x u) 2 S3 ( u £Re L ) ^ = O(Re^), and, since 
S3 d x u(x) /rfjf 1 , using the mesh size A from the relation (1), we arrive at the estimate 

1 . d n u(x) , A „ * 1 . „ . , A . „ , a n_1 „ i 
n\ ox n\ i] K n\ 

The relation (3) is essentially the basis for all numerical finite difference schemes used for the 
DNS of turbulence [9]. Indeed, we see that if a < 1, the first-order finite difference accurately 
represents the velocity derivatives. 

In spectral simulations of isotropic and homogeneous turbulence, one prescribes a suitable 
number of the Fourier modes to represent the velocity field. Usually, this number is chosen on 
the basis of the magnitude of the expected Kolmogorov scale tjk or the largest wavenumber 



kmax = 2k/vk- In the state-of-the-art simulations [10], [11], the cut-off is usually chosen such 
that k max = \/2N/3 on a grid of size iV 3 . 

In summary, the principal elements of Kolmogorov's phenomenology which have enabled 
these traditional estimates are the following: (a) the scaling exponents of the structure 
functions S n fl oc r^ n are given by the Kolmogorov values £ n = n/3; (b) the mean dissipation 
rate £ = u(diUj) 2 is constant and 0(1), as are the moments of the dissipation rate S n for all 
n; if the latter were not the case, one can define different Kolmogorov scales on the basis of 

n 

different moments of £; and (c) the "skewness" factors {d x u) n / (d x u) 2 2 = 0(1), independent 
of the Reynolds number; for, if this were not so, one can again define different Kolmogorov 
scales through odd moments of different order. 

The main point of the present paper is that there is a need to reexamine the traditional 
estimates in the light of modern developments in turbulent theory and experiment. We 
concentrate on isotropic and homogeneous turbulence but expect that the considerations 
hold for more general flows as well. 



2 Results for Intermittent Turbulence 

We are interested in the Navier-Stokes dynamics of incompressible fluids. In 1941, Kolo- 
mogorov derived the few exact relation of turbulence theory, presented here for an arbitrary 
space dimensionality d, as 

1 9 r d+1 S 3 , = (-1)^, 



r d+ig r v > d 

giving 5*3,0 = — d ^ 2 ) ^ r anc ^ ^s,o/^i,2 = 3. A dimensional generalization of this result, with- 
out however the analytical support, yields the Kolomogorov's (normal) scaling £ n = n/3. 
Recently [12], [13], some additional exact consequences of the Navier-Stokes equations have 
been derived. In combination with recent experimental results, we consider their conse- 
quences for intermittent turbulence. 

a. Dissipation scale as a random field We consider the moments of velocity difference 
(also called structure functions). Choosing the displacement vector r parallel to the "x-axis", 



we can define the structure functions S n ^ m (r) = (u(x + ri) — u(x)) n (v (x + ri) — t>(x)) n = 
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(S r u) m (S r v) n , where u and v are the components of velocity vector parallel and normal the 
x-axis, respectively. In the inertial range the velocity structure functions are Re- independent; 
that is, if the displacement r belongs to the interval n <C r <C L, then S ntm (r) do not involve 
any information about the dissipation scale. 

Modern experiments have revealed that Kolmogorov's result £ n = n/3 is almost cer- 
tainly incorrect and that £ n is a concave function of n — or the ratio £ n /n is a decreasing 
function of the moment number n. (See for example Refs. [14] for reviews and Ref. [15] 
for the most recent data.) Further, the form of structure functions is given by S 2n (r) = 
(u(x + r) — u(x)) 2n (2n— l)!!(eL)^(^) ?2n . The factor (2n— 1)!!, ensuring Gaussian statis- 
tics at the integral scale L, is a subject of a forthcoming paper, but it suffices here to say 
here that it has been recently verified in experiments and numerical simulations [16]. On 
the other hand, in the limit r — > 0, the analytic structure function is approximately equal to 
S2n(r) ~ (d x u(0)) 2n r 2n . Combining the two, we can define a natural dissipation scale of the 
2n th -order structure function [17]-[18] as 

In In /- [ 

V2n = {{d x u) 2n )^ n -^((2n- l)\\e- L-~ in Y^ n . (4) 

According to (4), the dissipation scales, which are expressed in terms of the moments of 
velocity derivatives, define a random field n. By a random field we mean here that the value 
of the length scale n depends on the order of the moment considered. It will be shown 
below that (4) is an approximation to a more accurate representation. Similar ideas were 
proposed earlier in Refs. [19]- [21] within the framework of multifractal theories. Writing 
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= [(2n — l)!!] 2n ~«2n ; and using the Stirling formula (n 3> 1), one obtains i 2n ~ (^) 3 for 
£n = n/3. This means that the effect of the factor (2n — 1)!! can be safely neglected. For 
anomalous exponents £ n < n/3, this factor is even closer to unity and does not modify the 
conclusions obtained below. 

b. Dissipation anomaly If the velocity field is different iable, we obtain S 3 (r) oc r 3 and 
drS^r) — > in contradiction with the Kolmogorov relation. This implies that the velocity 
field is singular in the limit of v — > and r — > (in that order), leading to the so-called 
dissipation anomaly. Here we first reproduce some details of Polyakov's derivation [22] of the 
dissipation anomaly for turbulence governed by Burgers equation and then outline similar 
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procedure for the Navier-Stokes equations. Consider the one-dimensional Burgers equation 



du du d 2 u 
dt dx dx 2 ' 



for which the energy balance reads as 



1 du 2 1 d o . . d 2 u 

2 dt 3dx v 'dx 2 

2' ^ u "" 1 ' 2dx± ~ ^dyi 



Introducing x± — x ± |, so that, | = ±^-, we can represent the energy balance equation 



as 



lim y ^ [ ±_ —+-—u(x+) 2 u(x-)+-—u(x-) 2 u(x + ) = I/ (_+_) u ( a:+ ) u (a;_))]. 

(6) 

We also have the identities: 

9. . NN o l r <9w(:r + ) 3 (9-u(a;„) 3 1 3 r 9w(x + ) 2 w(:r_) du(x_) 2 u(x . ), 

and 



where D, the dissipation contribution to the energy balance, is given by 



d 2 d 2 
D = v[u{x + )-^u{x+) + u{x-)-^u(x-)\. (9) 



Substituting these identities into the equation (6) and taking account of the fact that 
lim y ^,Q du ^ >3 = &U Q X J 3 , so that in the limit y — > all non-singular terms disappear by virtue 
of the energy equation (5), we are left with the balance between the singular (anomalous) 
contributions 

Id d 2 
l™ Qdy~( u ( x+ ^ ~ u ( x ~)) 3 = v [( u ( x +) ~ u ( x -))g^( u ( x +) ~ u ( x -))\- ( 10 ) 

This is Polyakov's expression for the dissipation anomaly derived for the Burgers equation 



[22]. Averaging (10) gives the exact relation {S y u) 3 = — V2£y where the dissipation rate 

du ~> 
dx> 
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We are interested in the Navier-Stokes dynamics of incompressible fluids, for which the 
energy balance equation (with the density p = 1) is written as 

1 du 2 1 _ , d 2 u 

and that for the scalar product u(x + |) • u(x — |) = u(+) ■ u(— ) can be written as 
du(+) ■ u(-) . . d . . . . , . 9 . . . . 

at + u(+) ' &^ u(+) ■ uH + u( -» ' &T uH ' u(+) = 

~W U - H - 1£r"' (+) + " uH ' a|- u(+) + u(+) ' sb uH| - 

It is clear that in the limit y — > 0, for which x± — > x, this equation gives the energy balance. 
Following Polyakov's procedure outlined above, let us consider the two identities: 



and 



<9 2 <9 2 
^(+)^2-^(-) +^(-)^2-^(+) = 

<9 2 <9 2 
-4(«i(+) - - «*(-)) + u <(+)^2-«i(+) + 



( 13 ) 

Similar identities for the pressure terms can be written easily. Substituting them into (11) 
and, as in the case of Burgers equation considered above, accounting for the energy balance, 
one has 

lim y ^o[~(u i (+) - u i (-))(u j (+) - Uj {-)f + i(-|_ Mi (+)^.(-) 2 + JL Ui (-) Uj (+) 2 ) = 

-4i/(u i (+) - Ui (-))^( Ui (+) - Ul (-)) + (Mtl _ ^1) . (u(+) _ U (-))].(14) 

This equation can be written in a compact form as 

d Id d 

lim y ^o[-—Su i \5 y u\ 2 + -{-—u i {+)uj{-) 2 + ^(-) M i(+) 2 ) = - 26 y u " M- 
y i H - ;^ 



where a = — Vp+z/V 2 u is the Lagrangian acceleration. The equation (14) is exact. Choosing 
the displacement vector along one of the coordinate axes and averaging (14), one obtains 



d d 2 4 

— 5u\5u\ 2 = 85ui—5ui = 2(5 y u i )d 2 (5 y u i ) = --£, 

where 5 y u = 5 y u -y/y. The pressure terms in and the second contribution to the left side of 
(14) disappeared by the averaging procedure. In general, we can choose a sphere of radius 
y « R — > around a point x and average (14) over this sphere. This causes the all 
scalar-velocity contributions to (14) disappear and the resulting equation can be perceived 
as a local form of the 4/3 Kolmogorov law. This fact has been realized before. Introducing 
the angular averaging, Robert and Duchon [23] and Eyink [24] locally expressed the relation 
(14) in terms of longitudinal and transverse velocity differences. We are interested in the 
order of magnitude estimates (see below), and restrict ourselves to (14). 

c. Relations between the moments In the isotropic and homogeneous turbulence, the 
Navier-Stokes equations lead to the following exact relations for structure functions. They 
were derived in [12] and [13] and experimentally investigated in some detail in Ref. [25]; see 
also Ref. [26]. The relations for different values of n are 

dS 2nfi d-l c _ (2n-l)(d-l) 



dr 



S 2n ,o = S 2n -2,2 + (2n - l)5 r a x (x)(5 r u) 2n - 2 . (15) 



Similar equations for all structure functions S njjn can easily be obtained from the equation 
for generating functions derived in [12]. 

d. The closure problem Equation (15), which includes both velocity and Lagrangian 
acceleration increments, is not closed and cannot be solved unless the relation between 
acceleration and velocity differences is established. It has been proposed in Ref. [17] that the 
local expression (14) written for the displacement magnitudes corresponding to the bottom 
of inertial range, i.e. in the limit y — > rj — > can be used as a closure. At the present time, 
this can be done only approximately. Since at the values of displacement y <C rj — > 0, the 
difference S y u ~ we can modify the lim operation in (14) as 

lim « lim , (16) 

y~ >0 y— >»y-^0 
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leading to the order-of-magnitude estimate 



.d(S v u) 3 d _. . . , 2 . d 2 . d5 v p(x) . . , 
hm y ^ v A — ^- h B—dyU{dyV) oc ud y u—d y u — — d y u ps d v ud v a x , (17) 



where A and 5 are undetermined constants. On extrapolating to the dissipation scale rj 
where all terms in the right side of (18) are of the same order, we derive the estimate [17] as 

v ~ i]5 v u = r)(u(x + i]) — u(x)). (18) 

The relation (18) tells us that each velocity fluctuation 5 v u is dissipated on its 'own' dissi- 
pation scale rj and the local value of the Reynolds number Rei = 0(1). This allows a simple 
physical interpretation that the dissipation processes at all levels n happen on "quasi-laminar 
structures" where the inertial and viscous terms are of the same order. In general, the higher 
the moment order, the more the intense events contribute, and the smaller the value of the 
corresponding dissipation scale. 

e. Dissipation scales and moments of derivatives The theory gives for the moments 
of Lagrangian acceleration a = — Vp + z/V 2 u the result that 

8„u (5 v u) 2 (S v u) 3 o Rc 
a x w -5- w w = (5 v uf -, (19) 

T v T) V U rms L 

where the turn-over time r T; m i]/5 v u. 

Below we will mainly discuss the equations for even-order structure functions, for which, 
if the displacement r is in the inertial range, the dissipation contribution to the increment 
of Lagrangian acceleration is negligibly small [17], [18]. For this case, we have 

0S 2nfi d-l Q (2n-l)(d-l 



+ S 2nfi = ± ^ J S 2n _ 2:2 - (2n - l)5 rPx {5 r u) 2n -\ (20) 

or r r 

where p x = d xP (x) and d denotes, as before, the space dimensionality. 

The relation (15) is valid for all magnitudes of displacement r L, including r — > 
77. Below, to simplify the notation, we will omit the subscript x in the re-component of 
acceleration a x . In this limit, treating (19) as a = lim r _,^((5 r -u) 3 /z/ and substituting it in (15) 
gives S2n r ^ S2n+ J ^ . On a scale r = r] 2n , writing S nfi oc A n rf™ , equation (15) gives 



r] n oc (21) 



For Kolmogorov turbulence with £ n = n/3 the formula (21) reads, as expected, as i] n = 

3 

r] K = LRe~4 which is n- independent. In intermittent turbulence, where the exponents can 
be well-described [12], [17] by the relation £ n fx 0.383 *n/(l + 0.05n), the relation (21) defines 
the Reynolds- number-dependent dissipation scales. As n — > oo, i] n — > LRe~ x . Thus, to 
resolve all fluctuations including the strongest, the computational work need to increase as 
i?e 4 , as already noted in Ref. [27]. In general, in the limit n — > oo, the relation (21) can be 
written as 



r] n fx LRe , 

so one may get a somewhat different estimate for the computational work than Re 4 , but 
the principal conclusion is inescapable that intermittency makes DNS more expensive than 
previously thought. 

Using the relations (18), (20) and (21), obtained by balancing various terms in the exact 
dynamic equations (14), (15), we can develop the multi-scaling algebra. For example, 

_ p„ p„ „.2 

FX (-^-fn S6n{r)6n) { _^_yn v l T „ (^Wn^a^ (22) 

with a 2n = 2n + ^ _|° n — ^ . With ^6 = 2 and £ 7 = 7/3, we recover Yaglom's result [28] 
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a 2 fx ^jf-. The intermittency corrections are readily found from (22). Recent experiments 
by Reynolds et al. [29] have lent strong support to this result. The formula (22) shows 
that the second moment of Lagrangian acceleration is expressed in terms of the sixth-order 
structure function evaluated on its dissipation scale i] e . To extract information about the 
fourth moment a 4 , we should have accurate data on S , i 2 (?7i2) which is very difficult to obtain 
in high-Reynolds-number flows. 

The moments of velocity derivatives are evaluated easily. In accordance with (18), we 
have 

„ ( ^)2n ~ ( (^0!)2n ~ ^ , (23) 

where d 2n = 2n H ^ 



It is important to stress that the first equality in (23) involves the averaging over two 
random fields u and r]. To perform this averaging, we have to either know the joint probability 
p(u,r),r) or use the functional relation between the fields given by (18). This leads to the 
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second equation in (23) and the final result. Since (d x u) 2 oc Re, the relation (3) leads to a 
new relation between exponents 

2^ = 6 + 1 

which agrees extremely well with experimental data. The relation (23) differs from proposals 
reviewed in Ref. [14]. 

f. The role of the fluctuations of the dissipation scale Let us reexamine the relation 
(4). In the limit r — > 0, the velocity field is analytic and can be expanded by Taylor series 



so that jf- m 5 r u/r. This gives (jr L r) 2n S 2n (r). When r — > rj — > 0, we have to evaluate the 



mean of the ratio (5 ri u/rj) 2n which is not a trivial task, since we are dealing here with the 
ratio of two random fields — unless the relation (18), which expresses the dissipation scale in 
terms of velocity field, is used. If, however, we incorrectly assume that the dissipation scale 
fluctuations are independent of those of the velocity field and neglect the step leading to the 
last equations in the right hand side of (23), it is possible to write the moments of velocity 
derivative as 

,8„u s 



( dxU )2n w ( ^_)2n w S 2n ( V2n ) / rj% oc Re^ , (24) 



where p 2n = g 2 _g 2 +"-r Equating expressions (23) and (24), we have 



^ 2 "~ 2 " = 2n+- ^ -, (25) 



&2n ~ &2n+l — 1 iin ~ 6n+l ~~ 1 

subject to the constraints £o = an d £3 = 1- The only solution to (25) is £ n = n/3. Since 
equation (25) is based on the first equality (23), which in general is incorrect, we can conclude 
that the source of anomalous scaling in hydrodynamic turbulence is the fluctuation of the 
dissipation scale field rj, which itself is strongly correlated the velocity field fluctuations via 
expression (18). This does not preclude a different situation from arising in other forms of 
turbulence, e.g., scalar turbulence generated by white- noise forcing [30]. 



It follows that (g ) 2 = lim^ %^ du ^ = - lim r ^ 2 £,u(x)u(x') oc (2 - fc)^" 2 . The 
higher-order derivatives are evaluated in a similar way to yield 



,d n u / d 2n £ 2 -2" e 2 -2 



n-l 



lim r ^ m \ —-S 2 {r)^r] 2 2 « Re 2 ^' 2 ) = Re^Re^- 2 . (26) 



dx n " mb '^" 2 \dr 2n 
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3 Implications for Numerical Methods 



2 « 



According to experimental data (see Refs. [25,15] for recent results), the exponent £ 
0.70 — 0.71 > 2/3 and as n — > oo, the terms in the expansion (2) for simulating the "typical" 
velocity derivatives can be estimated via 

(^) rms A"- 1 oc J Re^e^- 1 ), (27) 

with 7 = (— | — 1^2)) > 0- F° r £2 ~ 0.71, we find 7 ps 0.025. The accuracy of the numerical 
method in calculating the most intense velocity fluctuations can be estimated if, in the limit 
n — > 00, the expression 

i_ 

( | )2 „'"(A r . Kfe i fle ^ (28) 

is used instead of (d x u) rms . In the above equation, the mesh size A is defined by (1) and the 
expressions (23) for the moments of velocity derivative have been used. We see that when 
the Reynolds number is large, the high-order derivatives in the expression (2) dominate. 
This means that the DNS based on the mesh equal to the Kolmogorov scale becomes quite 
inaccurate. It is easy to check that accurate simulations of the largest fluctuations requires 
the resolution of the smallest scales which are 0(1 /Re). This means that the computational 
resolution scales as Re 3 and the computational work grows as Re 4 . 

In Refs. [19], it has argued that the intermittent nature of turbulence makes the size of 
the attractor smaller than the conventionally estimated, so the computational power needed 
becomes correspondingly smaller than the conventional estimate — not larger as just claimed. 
The rationale is roughly that the "interesting" parts of the flow occupy small volumes of 
space so any reasonable computational effort that focuses on those volumes is likely to be 
less expensive. This is also the spirit of adaptive meshing [31]. Even if the interesting parts of 
a turbulent flow are not space-filling, as discussed at length in Ref . [20] , we do not yet know 
how to track them efficiently in hydrodynamics turbulence. We also do not know if the part 
of the flow that contains the less interesting parts can be computed with greater economy. 
Nevertheless, it must be said that the present estimates apply to uniform meshing, which has 
been the most successful of the computing schemes until now. It should also be mentioned 
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that there is a specific suggestion [32] on the most singular structure in turbulence, which 
yields Re 3S , which is slightly different from Re 4 estimated in this paper. 



4 Dynamic Constraints on Sub-Grid Models for LES 

If the Reynolds number is large, the computational work involved in the numerical simulation 
of a flow is huge. It is interesting that at about the same time that DNS came into being, 
the idea of the Large Eddy Simulations (LES) was proposed by Deardorff [33]. The idea is 
very simple. Consider the Navier-Stokes equations 

d t u + UidiU = —Vp + ud 2 u; d, L Ui = 0. (29) 

We choose the mesh size A and define the so-called "sub-grid" velocity fluctuations ■u > (/c) ^ 
for k > 7r/A. The Fourier-transform of velocity field is defined as 

«(k) = u<(k) +« > (k), (30) 

so that 

«>(x) = / e ik ' x M > (k)rf 3 A;; u<(x) = / e^VCk^k. (31) 

The goal is to obtain the correct equation for the resolved scales u K {k) ^ in the interval 
< k < 7r/A. We decompose the field and write the equation for only the resolved scales as 

<9 t u< + uf ■ <9 iU < =SQ- Vp< + z/<9 2 u<, (32) 

where, for this particular formulation, the subgrid contribution is SQ = —uf ■ 9jU > — uf ■ 
(9jU < — uf ■ 9jU > . The LES equations are considered a success if the large-scale velocity 
fields (i.e., for k < 1/A) given by the Navier-Stokes equations (30) and by a model (33) are 
identical or close enough for all Reynolds numbers. 

There is, however, one problem. To derive the equation of motion containing only the 
resolved fields, one has to express all contributions to SQ, involving the sub-grid velocity 
fluctuations u > , in terms of u < , which is basically equivalent to solution of the proverbial 
"turbulence problem" . The model equation (33) is written in a generic form, but a similar 
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difficulty arises if, instead of the Fourier-space decomposition introduced above, the filtering 
or any other kind is used. 

The accurate LES model must satisfy the following dynamic constraints. The method 
developed in the Ref. [17] can be literally applied to the Navier-Stokes equations with an 
arbitrary right hand side and, defining the coarse-grained structure functions S< (r) = 



(S r u < ) n : we obtain, from (21), the result 
^ + ^^,o = (2W "^ (d " 1) ^-2,2 + (2" - l)(5 r (Sg x )-5 r p<)(5 r u<r-i. (33) 

The large-scale velocity fields obtained from DNS and LES can be identical S nj o(r) = S< (r) 
if and only if 



(5 r (SQ x ) - 6 rP <)(6 r u<y»- 2 = - WW 2 " -2 - (34) 

Similar constraints, coming from the equations for various structure functions S njTn can 
be readily obtained. It is impossible to demand equality of two random fields u and u< 
obtained from two different equations. The only criterion we can impose is that of statistical 
equality or, equivalently, constraint on all moments, namely S<(r) = S n (r). The relations 
(34), reflecting this necessary condition of the LES validity, must be satisfied. 

We wish to stress that these constraints are not dissimilar to S^f ~ S< m , often implied 
in the literature. Here S^^(r) are the structure functions evaluated from the velocity field 
obtained from LES. The velocity increment can be written as 5 r u = J u(k)e lkx (e lkr — 1), so 
that 

S 2 oc J E(k)(l - cos kr)dk. 

It is easy to see that if r << L, where L is the integral scale, and the energy spectrum 
decreases with k fast enough, the main contribution to the integral comes from the range 
where kr ~ 1. Thus the structure functions S n ${r-) probe structures on the scales of the 
order r and cannot differ strongly from the one obtained from the filtered field. 

Various model considerations, leading to expressions for SQ, have been suggested in the 
last forty years. Consider the example that follows from Kolmogorov's theory. If the role 
of the small scale fluctuations in the large-scale dynamics can be expressed in terms of 
effective viscosity use, then vsg ~ (eA 4 )^. Then, dropping the averaging sign (quite an 
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assumption!) and substituting a simple estimate coming from the energy balance, namely, 



e = vsgS^S^ = vscSfj, we derive the Smagorinsky formula [34] given by v S g = ay^S^A 2 , 
where a = 0(1). It is important that the resolved rate of strain is evaluated in terms of 
velocity differences on the computational mesh 



q < ( , 1 ufjx + ^-ufjx) ufjx + A^-ufjx) 

i>ijW - -( ^ + ^ ), (35) 



where i,j = 1, 2, 3. In this approximation, the Reynolds stress r^- = —UjU] m uSij ~ uscSfj. 
Equation (36) with the model for SQ defines a closed set of equations which can be used 
for LES. The analytically evaluated coefficient from Yakhot and Orszag [35] gives a ps 0.2, 
while the so-called dynamic method [36] gives something different. In all approaches, since 
the large-scale fields 5 r u < and 5 r u are statistically independent upon Reynolds number, the 
parameter a = O(Re ). Thus, this simple model is 

SQ « aA 2 V\S<\Vu < = 0(1). (36) 

Examining the relations (34) and (36), an interesting conclusion can be reached. If 
A <C r, one can assume statistical independence of all velocity differences <5 r -u < and 5a u k . 
Since SQ given by (34) and (35) depends on the velocity differences defined on the mesh size 
A as 



S r Sg(5 r u<) 2n - 2 w 5 r SG (5 r u<) 2n - 2 = 0, (37) 

we see that the Smagorinsky model satisfies the dynamic constraints, provided the pressure 
gradient differences in the filtered and unfiltered fields are close to each other. The validity 
of the dynamic Smagorinsky models in the range k << 1/A has been verified by large eddy 
simulations (A. Oberai, private communication 2005). However, as r — > A, S r SQ, S r p x and 
<5 r -u < are strongly correlated and, as a result, the model becomes invalid. This consideration 
is applicable to all low-order closures. 

This intrinsic failure of all existing LES models at scales comparable to the computa- 
tional mesh is well-known. At sufficiently low Reynolds numbers, LES give accurate results. 
However, with increase of Re the quality of the simulations deteriorates starting from the 
vicinity of the cut-off, propagating toward the larger scales. At this point one is forced to 
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increase the resolution. The reasons for this failure can be qualitatively understood as fol- 
lows. Consider LES at a relatively low Re on a fixed mesh A/Li = 71 where L\ is an integral 
scale of this particular simulation. Now increase the length scale of the flow L 2 » L 1: thus 
increasing the Reynolds number. If, in the first case, the number of the cascade steps for 
the energy flux to reach the mesh scale was say n 1; that in the second simulation is equal 
to Ji2 >> rii. Since the intermittency and deviation from the close-to-Gaussian statistics, 
experimentally observed at the integral scale, grows with the number of cascade steps, the 
contribution from the very strong velocity fluctuations at the "dissipation" scale A increases. 
As a result, the low order models that are successful in the close to Gaussian situations break 
down. In another scenario, let us increase the Reynolds number by increasing the mean ve- 
locity while keeping both the energy injection scale and the mesh size A constant. In this 
situation, the top of the "inertial" range will move into the range of scales which are larger 
than A, thus again invalidating the LES. 

A recent paper by Kang et al. [37] has demonstrated that, at the scales close to those of 
the mesh size, the probability density function p(S r u) computed from LES was quite close 
to a Gaussian while the experimental PDF showed broader tails, typical of intermittency. 
This means that the contributions from strong velocity fluctuations obtained from LES are 
underpredicted. Since the intermittent effects becomes stronger with increasing Reynolds 
number, we expect this difference to grow, thus invalidating the LES if the mesh size is 
also not modified. A very interesting example is given by the LES of the flow in a simple 
cavity reported by Larcheveque et al. [38]. It was shown that to correctly reproduce the 
experimental data on pressure fluctuations in a frequency range 100 < / < 2000Hz, the 
optimal cut-off of the large eddy simulations corresponded to Aj = lOOKHz. With decrease 
of A j, the quality of the results rapidly deteriorated. The present theory explains the failure 
of LES schemes with fixed mesh to describe the high Reynolds number flows as originating 
from the failure of low-order models in an all-important range r ~ A, this range being 
responsible for the energy cascade dissipation. At the present time, it is not clear how many 
constraints (35) must be satisfied to achieve accurate LES, but we believe that the number 
must grow with the Reynolds number. 
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5 Conclusions 



For many years, intermittency and anomalous scaling in three-dimensional turbulence were 
considered major challenges for theorists. Recent developments of the multifractal theory 
and its dynamic formulation led to description of intermittency in terms of an infinite number 
of dissipation scales (ultraviolet cut-offs). It was shown that strong velocity fluctuations are 
dissipated on scales that are much smaller than that estimated from Kolomogorov's theory. 
In this paper, we have attempted to make a connection between the theory of anomalous 
scaling and numerical methods. 

One conclusion that follows from this connection is that to simulate all fluctuations, 
including the strongest ones, the computational demands scale as i?e 4 , and not as Re 3 as 
traditionally deduced according to the Kolmogorov theory. To achieve the full DNS of tur- 
bulence, including the strongest small-scale velocity fluctuations, one has to use resolutions 
high enough to produce an analytic interval of structure functions, where S n ~ d x u(0)) n r n . 
Analyzing the results of various numerical state-of-the-art DNS, we have discovered that this 
criterion is satisfied only for n < 4. This is not sufficient to accurately simulate the velocity 
derivatives. 

A second comment concerns the Large Eddy Simulations. An infinite number of dynamic 
constraints on a correct subgrid model has been derived from the exact relations for structure 
functions. Due to the Galilean invariance, the subgrid scales cannot influence the advective 
term in the Navier-Stokes equations, provided the subgrid scale A/r — > 0. However, it is 
clear from analyzing the equations of Section 4 that the subgrid model cannot be reduced 
to a low-order viscosity expression, but must include high-order nonlinear contributions that 
do not vanish at the scales close to the mesh size. 

Thus, while accurate DNS are possible if the resolution requirements are met and powerful 
enough computers are available. However, due to the basic theoretical problems, derivation 
of an accurate and theoretically justified subgrid model, valid at very high Reynolds numbers, 
remains a major challenge. 

It is worth pointing out that we have considered homogeneous and isotropic turbulence. 
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The situation with wall flows is even more complex. There, turbulence is mainly produced in 
the vicinity of the wall where acceleration and turbulence production are highly intermittent. 
Recent DNS by Lee et al. [39] have demonstrated strong intermittency and the Reynolds 
number dependence of the few first moments of Lagrangian acceleration near the wall, sharply 
peaking at the reduced normalized distance y + ~ 2.5. At present, we do not know how to 
model this near-wall phenomenon that is largely responsible for turbulence production. 

We wish to conclude on a "positive" note. The fact that the structure functions S271 ~ 
(2n — l)!!(-jr)£ 2n means that the velocity distribution is close to the Gaussian near r = L, 
and the intermittency is weak or nonexistent. It follows that simple, semi-qualitative re- 
summations of the expansions in powers of the dimensionless rate-of-strain are much less 
problematic there. Thus, the derivation of the VLES or time-dependent RANS appears to 
have a brighter future. 
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